Self-Organisation in 2D Swarms. 
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We undertake a systematic numerical exploration ol self- organised states in a deterministic model 
of interacting, self-propelled particles in 2D. In the process, we identify new types of collective 
motion, namely disordered swarms, rings and droplets. We construct a "phase diagram", which 
summarizes our results as it delineates phase transitions (all discontinuous) between disordered 
swarms and vortical flocks on one hand, and bound vortical flocks and expanding formations on 
the other. One of transition lines is found to have a close analogy in the temperature-driven gas- 
liquid transition in finite clusters with the same interparticle potential. Furthermore, we report on 
a novel type of flocking which takes place in the presence of a uniform external driver. Altogether, 
our results set a rather firm stage for experimental refinement and/or falsification of this class of 
models. 
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Populations of self-propelled organisms tend to orga- 
nize in remarkable aggregate formations Schools 
of fish, flocks of birds, swarms of bees or locust are 
among the more familiar examples. Less familiar per- 
haps, though equally prevalent, are the organized states 
achieved by microscopic organisms: amoeba aggregates 
Q , bacterial colonies 0, 0] , and swarms of Daphnia @ . 
Such striking self-organized collective phenomena stimu- 
lated natural extensions of trusted models and tools of 
equilibrium statistical mechanics to systems of coupled 
self propelled particles. Proposed models have varied 
in degrees of idealization and complexity of the descrip- 
tion of swarms, their environments and the interactions 
between them. They can be broadly divided into con- 
tinuous [sHp and discrete self-propelled particle (SPP) 
models 10l4l3|. both coming in deterministic as well as 
stochastic varieties. In these pioneering studies, mod- 
els were mainly analyzed for their ability to display cer- 
tain observed phenomena (transition from swarm to flock 
H)[I3|> vortical swarms 12], as well as the potential 
universality of certain features (scaling behavior 11| , in- 
dividual distance (H)]). More often than not, they were 
not pushed to predictive risks that would at best qualify 
and refine them, at worst falsify them. While these ex- 
ercises are crucial to identify generic sufficient conditions 
for the occurrence of observed states, they leave open the 
question as to whether or not these conditions necessarily 
obtain in the living system under investigation. In cases 
where such systems are open to detailed experimental 
investigation Q, it proves useful to explore the richness 
of behavior sustained by classes of models, with a view 
to proposing tracks for spatio-temporal evolution that 
might result under (slow) changes in the experimental 
setting. This is precisely the objective we had in mind as 
we embarked on our explorations, and shall dedicate the 
rest of this letter to describing the working model (dis- 



crete, deterministic, SPP model), summarizing our ex- 
tensive numerical experiments, and highlighting our sig- 
nificant, and eminently testable, results: a- novel types of 
self-organisation namely the disk, ring, polarized vortex, 
droplet and expansion states; b- a phase diagram which 
distills, on a reduced parameter plane, the various states, 
and the transitions between them; c- an intriguing route 
to flocking, which obtains in the presence of an external, 
uniform, force field. 

Given the inherent complexity of our undertaking, we 
opted for a simple, but fairly versatile, deterministic SPP 
model, flavors of which have been examined by Levine et. 
al. [13, in their quest for vortices, and Edelstein-Keshet 
et.al [151 ] in a study of individual distances in swarms. 
The particles in this model are identical in their mass, 
and in the nature of the forces they feel and generate. 
Their self-propulsion is mimicked with an acceleration 
of constant magnitude, acting along the direction of mo- 
tion. They are coupled via a double exponential potential 
force field, which attracts at large distances, and repels 
at small distances. Furthermore, they are subjected to 
a drag force, which captures the reaction of the viscous 
medium in which they move (and which is here assumed 
linear in the velocity). Newton's equation for the zth 
particle reads 



dvi „ 
at 



(1) 



where m, f\ and vl are the mass, position and velocity of 
the ith particle, a the magnitude of the self-propulsion 
force, which acts along the direction of motion Vi, and 7 
the friction coefficient. The potential energy <p is given 
by 

N 

(j)= ^2 W r exp(-\f' i ,-f}\/l r )-W a exp(-\r i -f' j \/l a ) (2) 

i,i=i 
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where N is the population size, and W r , W a , l r and l a 
determine the strength and range of repulsive and at- 
tractive forces respectively. It proves useful to work with 
the dimensionless variables r' = r/l r , v' = v/v t , and 
t' = t/r, with vt = — (the "terminal velocity"), and 
t = l r /v t = l r j/cr. In the primed variables, the equation 
of motion becomes 



R 



N 



3 = 1 



(3) 

with = [fj — r^/lfj — r^|. We recover four di- 
mensionless parameters: A = l r /l a (repulsive over at- 
tractive length scale, < 1), P = W a /crl a (attraction 
over self-propulsion) , Q = W r /al r (repulsion over self- 
propulsion), and R = ma/(l r j 2 ) (a Reynolds- like num- 
ber for the flow) , which together with the population size 
N, fully specify a swarm in our model. 

In what follows, we are primarily concerned with the 
behavior of 2D swarms. Initially (and except when oth- 
erwise specified) , particles are uniformly distributed in a 
square (its size being chosen so that the initial swarm is 
bounded [HI]), with randomly oriented initial velocities, 
and speeds uniformly distributed between and the ter- 
minal velocity. Populations ranging in size between 100 
and 1000 members were considered. With these initial 
conditions, and the help of a 4th order, adaptive, Runge- 
kutta scheme |17( , we thoroughly explored the parameter 
space. Simulations were allowed to run till the system's 
energy (and wherever relevant, mean velocity and den- 
sity) relaxed to near steady configurations, on which our 
classification is based. Results are viewed in the center 
of mass frame, the origin of which, though interesting in 
its wanderings, is not a serious concern of this inquiry. 

Broadly speaking, these random initial conditions con- 
verge, to a center manifold (think of it as the terminal 
velocity manifold), on which the swarm relaxes to dis- 
tinct regimes, which though not exhaustive of the rich 
behavior allowed by this model, constitute the skele- 
ton on which a complete description will be eventually 
fleshed out. Before displaying, then analyzing, a typical 
(two dimensional) space in which these regimes find there 
proper home, we survey their main qualitative properties. 
We draw our results from experiments with 1000 parti- 
cles, while noting that the existence and main qualitative 
properties of these regimes are independent of the size of 
the population. 

We start with the ring state, shown in Fig. QJa), a 
highly regimented configuration, in which particles ro- 
tate and counter-rotate on a perfect circle. Next in line, 
is the polarized vortex, shown in Fig.fljb), in which nearly 
all particles travel in the same direction 18| . A hith- 
erto unsuspected formation is the droplet state, Fig. QJc- 
d), in which the swarm breaks up into a necklace of 
drop- like flocks. Droplets rotate clockwise and counter- 
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FIG. 1: Steady states, for N = 1000, A = 2/3, and P = 0.1: 
(a) Ring State, Q — 0.1 and R = 3. (b) Polarized vortex, 
Q = 0.21 and R = 0.0975. (c) and (d) Droplet states, with 
Q = 0.1575, R = 3.75 and Q = 0.1125, R = 7.5 respectively, 
(e) Vortex state, Q — 0.15 and R = 5.625. (f) Disordered 
Disk State, Q = 0.1875 and R = 0.0375; (g) Expansion State, 
Q = 0.2625 and R = 0.375. 



clockwise, with near constant angular velocity, hold- 
ing tight through repeated mergers with fellow droplets 
traveling in the opposite direction. Bona fide vortices, 
on which particles are nearly split even between pro- 
grade and retrograde circulation, are shown in Fig. QJe) . 
Fig. [Iff ) shows an example of the disordered disk state, 
in which particles move chaotically, while self-confined 
within a disk like region. Lastly, the expansion state, in 
which self-propulsion dominates over attractive forces to 
promote the outward explosion of the swarm, is shown in 
Fig. Eg). 

Perhaps more important than the identification of 
novel self-organized states, is the realization that the 
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FIG. 2: The phase diagram, for N = 1000, and A = 2/3. 
P = R/NP is plotted against a = Q/P. The dotted a = 1 
line, the line of rings, separates the multi-swarm (attraction 
dominated) region to its left, from the swarm dynamics re- 
gions to its right. Circles mark transitions from disordered 
disk to polarized vortex state; triangles mark transitions from 
ordered to expansion state. 

five dimensional parameter space (the four dimension- 
less model parameters, plus N) reduces, once the dy- 
namics settles on the terminal velocity manifold, to a 
three dimensional space, in which the variety of re- 
laxed, self-organized configurations, is itself organized 
in a phase (morphology) diagram of sorts. This space 
is spanned by A = l r /l a , a — Q/P — W r /(XW a ), and 
(3 = R/NP = mVf/(NXW a )- We verified, with extensive 
numerical explorations, that models with the same (A, 
a, /3) relax to qualitatively similar configurations, that a 
steady configuration in the (A, a, j3) plane unfolds into a 
hyper-surface of similar configurations in the five dimen- 
sional parameter space of the model. 

It was noted [19j that the Morse potential (Eq. [5]) with 
parameter values resulting in self-organized states is not 
i/-stable in the language of statistical mechanics [2fjj |. 
-ff-stability is related to the existence of a well-defined 
thermodynamic limit with finite intensive local densities. 
The system which is not -ff-stable is called catastrophic. 
Indeed, the necessary criterion to ensure this stability is 
that the integral over space of the pair interaction poten- 
tial be non-negative. Since 

J c/)(r)d 2 r = 2-kWJI ("A 3 - l ) 

the criterion is not satisfied and the system is catas- 
trophic for a < A~ 3 . In the catastrophic regime the 
average potential energy per particle in large but finite 
systems is proportional to N. 

The reduced set of parameters has, therefore, simple 
physical meaning: a and A together define the shape 
of the potential, and, in particular, the degree to which 



the system is away from .ff-stability. As a increases to- 
wards the value of A -3 the system approaches the nor- 
mal, non-catastrophic regime. On the other hand, j3 has 
the meaning of the ratio of the typical kinetic energy to 
the potential energy per particle (which is dominated by 
attraction and is of order NW a ). 

In this reduced parameter space, one could survey 
steady states, and neighbouring configurations, by vary- 
ing three of the model parameters (say W a , l a and 7 
or a) and keeping the others (W r , l r , N and a or 7) 
fixed. Thousands of computing hours went into clarify- 
ing the salient features of the resulting phase diagram, 
the invariant structures of the model. These features are 
clearly delineated on constant A slices, where the relaxed 
dynamics typically splits into four regions, as apparent 
in the particular instance in Fig. [5] a) a region occupied 
by the disordered disk states, for smallish j3 (attraction 
and viscosity dominated models); b) a region occupied 
by unbound, expanding states, for largish /? (propulsion 
dominated models); c) an intermediate region where all 
organized states live (i.e. rings, droplets, polarized vor- 
tices and regular vortices). This region is bounded be- 
low by the critical phase transition line (CPT1) between 
disordered and ordered states, and above by the criti- 
cal phase transition line (CPT2) between ordered and 
expanding states. It is bounded to the left by the line 
a = 1, or W a /l a = Wb/lb, a line of ring states, the radius 
of which increases with increasing j3, all else being held 
fixed, till the point (past CPT2) where the dominant self- 
propulsion fragments the ring into an expanding state; d) 
a region of irregular, multi-cluster, and unsteady configu- 
rations, which is obtained in attraction dominated mod- 
els (to the left of the ring line, and which we keep out of 
our humble focus on coherent, relaxed swarms). 

CPT1 and CPT2 meet at a crossroads between the 
disordered disk, the polarized vortex and the expanding 
state; CPT2 stretches beyond that meeting point into 
regimes which are held together by viscosity's countering 
of the mutual repulsion of particles. For a > 1, and /3 just 
above CPT1, polarized vortices are always observed. In 
fact, just above CPT1, all particles are circulating in the 
same direction. Increasing j3, at constant a, the number 
of counter rotating particles grows steadily till the po- 
larized vortex is destabilized, in three possible ways, de- 
pending on the value of a: 1) The polarized vortex breaks 
up into droplets (mostly two) rotating in the same direc- 
tion, with a small fraction of particles moving erratically 
in the opposite direction. Increasing /3, the number of 
droplets increases, with nearly as many droplets rotating 
in one direction as the other. Looking at larger values 
of j3, the droplets can either go directly to the expansion 
state or, and for a small a interval, they become more 
elongated and pass through the vortex state before break- 
ing up and expanding. 2) The polarized vortex morphs 
directly into a regular vortex and then to an expand- 
ing state. 3) The polarized vortex transitions directly to 
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the expansion state; here, and very close to CPT2, we 
find mixed states-that is states in which a small frac- 
tion of particles expands while the other forms a special 
vortex 2l| . Transitions among the organized states are 
harder to identify numerically, better handled through a 
careful stability analysis in a mean field approach to the 
problem, and are only roughly determined in the cur- 
rent study [Hf. On the other hand, transitions across 
CPT1 and CPT2 are discontinuous, and reminiscent of 
the first order phase transitions of equilibrium statistical 
mechanics. To better clarify this character, we define an 
order parameter L =< r\9\ >, the bracket indicating an 
average over all the particles in a given steady state. L 
is nearly equal to zero in the expanding state, tends to 
— « 0.636 for a completely random state, and 1 for vor- 
tex like dynamics. Following L with /3, for fixed a, we 
observe the system in Fig. [3] undergoing the two sharp 
transitions, from random to organized, and from orga- 
nized to expanding. For a > 2.7, the swarm transitions 
directly from the disordered to the expansion state. 
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FIG. 3: The transition among CPT1 and CPT2 are shown. 
The order parameter < r\8\ > is plotted against f3. The graph 
was obtained for N = 1000, A = 2/3, and a = 1.275. 

To clarify the nature of some of the transitions ob- 
served, we undertake a study of a thermodynamic anal- 
ogy of the self-propelling system. 

It is clear from the phase diagram that one of the con- 
trol parameters, /3, has the meaning of the ratio of the 
typical kinetic energy to the potential energy per parti- 
cle. In the thermodynamic analogy, we consider a system 
with exactly the same interaction, but the kinetic energy 
is controlled by temperature instead of the self-propulsion 
to viscosity ratio. Consider the canonical ensemble de- 
scription of a system of N particles with pair interaction 
specified by Eq. [5] enclosed in a circular box of radius 
R and kept at temperature T. Thermodynamic proper- 
ties are treated in the self-consistent field approximation. 
According to the standard formulation of the variational 
mean field theory [23| the density is related to the mean- 
field potential by 



where A is determined by normalization 

N = n(f)d 2 r 



and the mean- field potential is given by 
fmf(r) = / n(r')(j)(f' - f)d 2 r' 



Axial symmetry was assumed, and self-consistent field 
equations were solved numerically by an iterative pro- 
cedure for various temperatures and interaction param- 
eter a at fixed A=2/3. In a fairly broad temperature 
range two stable solutions exist: one with an almost uni- 
form low density, corresponding to the gas phase, and the 
other with a large density near the center and extremely 
low density at the periphery, corresponding to the liquid 
droplet phase. Intriguingly, stable droplet solutions exist 
only within the catastrophic regime a < (3/2) 3 . 

In order to analyze thermodynamic stability we calcu- 
late the Helmholtz free energy for each solution (phase) 
in the same self-consistent field approximation: 

F = |J n(r)n{r')4>(\r - r l \)d 2 rd 2 r'+ 

+T J n(r)\n (n(r)/N)d 2 r 

The velocity-related contribution is the same for both 
phases and will not be considered explicitly. The solu- 
tion with the lower free energy corresponds to thermody- 
namic equilibrium while the other solution can be associ- 
ated with a metastable phase. Upon approaching certain 
temperature, a solution may lose it iterative stability and 
eventually disappear. In the language of thermodynam- 
ics, this would mean reaching the spinodal line. Temper- 
ature dependence of the two branches of the free energy 
is shown in Fig.[4|a). 

The fact that two branches of the free energy cross is 
typical of classical first order transitions. However, in 
contrast to standard situations simultaneous coexistence 
of two phases is impossible. Each of the two solutions de- 
scribes the system as a whole while, say, a 50-50 mixture 
of uniform gas and liquid droplet is not a self-consistent 
solution at all. This is why the transition temperature is 
defined from the condition that the two phases have equal 
total free energies rather than from equality of chemical 
potentials. We recall that another distinction of the cur- 
rent situation is that a standard thermodynamic limit is 
not very meaningful. However, from a pragmatic point 
of view, first order transitions in finite systems can still 
be defined if there are two distinct states with the free 
energies that are large on the T scale and cross as func- 
tions of some control parameter. Phase transitions in 
an isolated macromolecule are conceptually close to the 
situation studied here, 



n(f) = Aexp 



T 



24] 



We make a connection with the dynamic system by 
equating mv 2 and T (Boltzmann constant is 1) since 
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FIG. 4: (a) Temperature dependence of the two branches 
of the free energy: liquid droplet phase (solid line) and gas 
phase (dotted line) for N = 1000, A = 2/3, a = 2.85. The 
lines terminate at respective spinodal points. T is expressed in 
W a units. (b)Phase diagram of the thermodynamic system 
for TV = 1000, A = 2/3, box radius R = 30Z a , in the same 
coordinates as in Fig [5] The equilibrium transition line is 
solid, spinodal lines are dashed. The liquid droplet phase is 
stable below the transition line, and metastable in the range 
between the solid and the upper dashed lines. 



the system is 2-dimensional, so that (3 — T/NXW a ■ 
The transition line in the (a, Log/3) coordinates is pre- 
sented in Fig. IDJb) together with two spinodal lines for 
iV=1000, A=2/3, and R = 30l a - Separate calculations 
with a different number of particles iV=500 confirm that 
our choice of j3 and a as scaling parameters yield a uni- 
versal A-independent phase diagram. There is a weak 
(logarithmic) effect of the box radius R on the transition 
temperature which is naturally related to the entropy of 
the gas phase, but further elaboration of the thermody- 
namic model is the subject of future work. 

It is clear that the gas-droplet transition line resembles 
closely the CPT2 line separating the expanding and the 
ordered states. This gives a strong support to the claim 
that the transitions in the dynamic system are indeed a 
direct analogy to phase transitions. Analysis of the dy- 
namic analogy of the metastable thermodynamic states 
would involve a study of the attractor basins in the phase 
space and is beyond the scope of the present paper. The 
thermodynamic system does not support dynamically or- 
ganized states, nor do we observe a transition analogous 
to the CPT1 line. A system of self-propelled particles 
can be highly organized in the momentum space which is 
impossible for a classical system of interacting particles 
in the canonical ensemble where Maxwell's distribution 
is inavoidable. Note, however, that the current state- 
ments concerning the thermodynamic system are based 
only on the results of the self-consistent field approxima- 
tion, while the true ordering in the coordinate space may 
be more intricate. 

We conclude our discussion of critical phase transitions 
with a brief mention of an interesting phenomenon, which 
is analogous to the notorious swarm to flock transition, 
and which, unlike the self-promoted transitions seen in 
this and other instances in the literature, occurs in the 
presence of an external, uniform, force field. Fig. [5] (b) 
shows a polarized flock in which all particles are traveling 
with the same velocity, a state which resulted by subject- 
ing a disordered disk state of our model (seen in Fig. [S] 



(a)) to a constant, external force field. Such a polarized 
flock does not emerge for any field strength. In fact, a 
discontinuous transition, from random to polarized, oc- 
curs as the strength of the external field, f ext , is increased 
past a critical value (which naturally depends on the ini- 
tial phase of the swarm (i.e. on the swarm's A, a, and 
0). The center of mass speed: u) = \^2f =1 vl\/ (Nv max ) 
(normalized by v max , the maximum speed in the swarm), 
is a suitable order parameter for the transition from dis- 
ordered disk (to — 0), to polarized flock (uj — 1), which is 
clearly illustrated in Fig. [5jc) . 



FIG. 5: The transition from disordered swarm to an aligned 
flock: a) Shows the initial disordered state; b) Shows the 
steady solution after applying the external field (fext = 
0.00115); c) shows the 1st order phase transition where we 
plot the order parameter u vs the external field strength (nor- 
malized by self-propulsion a). The graphs were obtained for 
N = 1000, A = 2/3, a = 1.875 and f3 = 5.25 x 10" 5 . 

In summary, we have succeeded in constructing a phase 
diagram that captures, in one fell swoop, the global struc- 
ture of steady 2D configurations, in a deterministic model 
of coupled self-propelled particles. Through a systematic 
numerical exploration of the model's parameter space, 
we have identified novel types of collective behavior such 
as rings, droplets, special vortices, disordered disks, and 
expanding states; we have shown that transitions from 
disordered to ordered, then from ordered to expanding 
state, have the structure of discontinuous, "first order", 
phase transitions (in line with the results of [25j |. which 
attribute the continuous transitions reported by to 
finite-size effects); we have supported the phase tran- 
sition analysis by studying a thermodynamic analogy 
where a temperature-driven transition from a gas to a 
condensed liquid droplet is observed and matches closely 
the transition from ordered to expanded state; finally, we 
have shown that a disordered swarm, when subjected to 
an external uniform field, transitions to a polarized flock, 
once again in a discontinuous manner. The genericity of 
discontinuous phase transitions, the stability and evolu- 
tion of organized states, the surface tension and shape 
dynamics of swarms, are the subject of ongoing work. In 
the mean time, the constructed phase diagram offers a 
promising geography for experimental verification, and 
refinement, or at worst, falsification of this class of mod- 
els. In particular, the same bacteria which stimulated an 
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earlier investigation of model vortices 12[ , could perhaps 
be stimulated into transitions from disordered motion, to 
polarized vortex, to droplets, or vortex, by judicious se- 
lection and control of population size, medium resistance, 
concentrations of chemo-attractants/repellents, external 
drivers... etc. As evident in Fig. [51 an order of magnitude 
increase in population size (decrease in /3) is expected to 
bring about the disintegration of a colony, initially orga- 
nized in a polarized vortex, into a disordered disk, and 
this over a range of interaction potentials (of as); simi- 
lar effects could result from a slowing down of particles 
in the colony (through increased friction and/or reduced 
self-propulsion); an increasing concentration of chemo- 
attractants (a decreasing a) may pus h an initially vorti- 
cal colony closer to a ring state 26]. We look forward 
to exchanges with experimentalists around this model's 
phase diagram, in the hope that similar systematic (not 
to say exhaustive) explorations of macroscopic phases of 
collective motion in the lab, together with the model im- 
provements that they will surely stimulate, may eventu- 
ally pave the way to a refined characterization of the no- 
toriously elusive (at times microscopic) mechanical prop- 
erties of coupled, self-propelled living organisms. 
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